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Abstract. We describe an assembly of numerical tools to model the output data of the Planck satellite <Taubedl20 04'). These 
start with the generation of a CMB sky in a chosen cosmology, add in various foreground sources, convolve the sky signal with 
arbitrary, even non-symmetric and polarised beam patterns, derive the time ordered data streams measured by the detectors de- 
pending on the chosen satellite-scanning strategy, and include noise signals for the individual detectors and electronic systems. 
The simulation products are needed to develop, verify, optimise, and characterise the accuracy and performance of all data 
processing and scientific analysis steps of the Planck mission, including data handling, data integrity checking, calibration, map 
making, physical component separation, and power spectrum estimation. In addition, the simulations allow detailed studies of 
the impact of many stochastic and systematic effects on the scientific results. The efficient implementation of the simulation 
allows the build-up of extended statistics of signal variances and co-variances. 

Although being developed specifically for the Planck mission, it is expected that the employed framework as well as most of 
the simulation tools will be of use for other experiments and CMB-related science in general. 

Key words, cosmic microwave background - large-scale structure of the Universe - methods: numerical 



1. Introduction 

The simulation of realistic data streams is a very impor- 
tant task in the preparation for the Planck satellite mission. 
Scientifically, simulated data are required for a realistic assess- 
ment of the mission goals, and for planning how these goals 
can best be achieved. The performance of data-analysis algo- 
rithms and software needs to be tested against the scientific 
requirements of the mission. This pertains to methods for re- 
moving systematic noise from the maps, separating the multi- 
frequency data into different physical emission processes, de- 
riving power spectra and higher-order statistical measures for 
CMB temperature and polarisation fluctuations, identification 
of "point" sources like distant galaxies and galaxy clusters, and 
many more. In order to develop, test, and calibrate methods 
for pursuing these scientific goals, full-sky maps need to be 
produced containing the CMB and all known components of 
foreground emission in the nine Planck frequency channels be- 
tween 30 and 857 GHz, at the resolution of 5' and the sensitiv- 
ity of AT/T ~ 10"^ to be achieved by the Planck instruments. 
These maps need to be "observed" following the scanning law 
foreseen for Planck and using realistic beam shapes, filtered 
with the frequency and time response of Planck's detectors, 
and deteriorated by random and systematic instrumental noise 
before realistic scientific assessments can be made. 
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Also important is the usage of simulated data for the devel- 
opment of data-analysis techniques for the mission. At a resolu- 
tion of 5', and assuming sufficient oversampling of the beams, 
the full sky has on the order of 5 x 10^ pixels. Full-sky maps 
in all frequency bands and physical components (like fore- 
grounds) need to be efficiently handled and stored. There will 
be 74 detectors on-board Planck, 52 on the High-Frequency 
Instrument (HFI) and 22 on the Low-Frequency Instrument 
(LFI). During one year of the mission, these detectors will pro- 
duce on the order of 10" measurements along the slowly pre- 
cessing, stepwise-circular scanning path, from which the fre- 
quency maps need to be produced. Storing the data efficiently, 
preparing them for fast access, and handling them within the 
data-analysis software pipelines needs to be thoroughly tested 
using data simulated at the correct temporal and spatial sam- 
pling rate. 

Simulated data are also important for assessing mission op- 
erations and the technical planning. Depending on the scanning 
law, i.e. the path along which the satellite moves and scans the 
sky, the noise characteristics of the maps will be more or less 
anisotropic because some regions on the sky will be observed 
much more often than others. At a fixed scanning law, the posi- 
tion of the detectors in the focal plane, their orientation relative 
to the optical axis, and the orientation of the optical axis rela- 
tive to the spin axis of the satellite further affect the scientific 
goals by changing the noise characteristics, the coverage pat- 
tern on the sky, the quality of polarisation measurements and 
more. 
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Finally, it is an important side-aspect of simulated data that 
they can be used to investigate the advantages and disadvan- 
tages of different models for data organisation and storage. 
With single data sets exceeding the size of multiple GBytes, ef- 
ficient data handling and processing requires that these data be 
stored in highly sophisticated ways, e.g. in data bases such that 
sequential and random access to parts of the data is possible. 
As a European-based mission with many participating coun- 
tries, Planck data analysis will operate in a distributed fashion 
(i.e. in several geographically separate data processing centres), 
which places further stricter requirements on the way data are 
archived and retrieved. 

The arguments above demonstrate that the simulation 
pipeline is essential for, and an ir itegral part of, both the HFI 
and LFI Data Processing Centres jPasian & Svgnei2002h . 

Even before the announcement of opportunity for the 
Planck instruments was addressed, and actually during the 
phase-A studies for the mission, several simulations were made 
to understand the possible science outcome from Planck and to 
define the instrumental setup. The need was felt to put on the 
same footing such simulation effort, and to extend within a co- 
ordinated environment. 

These considerations led to the construction of the Planck 
simulation pipeline which began quite early in the project, 
approximately six years ago. At the time there were numer- 
ous simulation activities within many groups involved in the 
project, but those simulations were incomplete, aiming at in- 
vestigating a multitude of specific and different goals. They 
were widely used by different groups throughout the Planck 
consortia, but produced data in non-standardised and non- 
portable formats, and were not optimised with respect to their 
consumption of computer resources. 

Starting with such existing simulation programs, the Planck 
simulation pipeline was built according to the following design 
guidelines: 

- The pipeline must be modular, i.e. different simulation 
tasks must be carried out by separate programs, called mod- 
ules. This simplifies extending the pipeline by adding fur- 
ther simulation tasks as they become necessary, and it also 
facilitates testing and maintaining the pipeline and its mod- 
ules. 

- Data must be exchanged between modules in a sim- 
ple, standardised, and machine-independent way. For the 
pipeline to be modular, the input-output interfaces of the 
modules are crucially important. For the interfaces to be 
simple and platform-independent, it was decided that data 
be exchanged in the FITS ' format. It is foreseen that most 
of the data transfer between Planck simulation and data- 
analysis modules will ultimately be done through a data 
base, but this does not change the principle of exchanging 
data through simple and platform-independent interfaces. 

- Modules can be written in any suitable programming lan- 
guage. Since the only contact of a module to its environ- 
ment is established through its interfaces, the internal im- 
plementation of the module is of no importance. Instead, 

' |http://archive.stsci.edu/fits/fits_standard/| 



allowing programmers to choose the language they found 
most suitable or convenient allowed the pipeline to be 
rapidly extended. 

Following these criteria, the construction of the Planck sim- 
ulation pipeline proceeded in four steps. 

- First, simulation programs were collected within the Planck 
consortia, converted to modules in the sense that they 
were enclosed within suitable interfaces, and augmented 
by missing modules specifically written to complete the 
pipeline. In that way, it was possible within three months to 
produce the first simulated time-ordered data streams from 
full-sky temperature maps of the CMB and several fore- 
ground components, which were convolved with realistic 
beams along the scanning path, filtered according to the fre- 
quency response and appropriately time-sampled, and had 
random noise added. 

- Second, the function of the pipeline was extended by 
adding standard-compliant modules necessary to simu- 
late the observation of polarised signal components, point 
sources, systematic noise, extended beam wings, and the 
CMB dipole, taking the exact motion of the spacecraft into 
account. Massive production of simulated data was started 
(e.g. time-ordered data were produced for different detec- 
tors, scanning strategies and combinations of foregrounds). 

- In a third phase, the Planck pipeline was consolidated by 
improving all individual modules in terms of their algo- 
rithms and their implementation, and adding modules sim- 
ulating higher-order effects like the emission of moving 
point sources or the detailed motion and rotation of the 
spacecraft. In parallel, a large part of the code base of the 
pipeline was ported from several programming languages 
to C+-I-. 

- Finally, the entire code was optimised in the sense that 
memory consumption was drastically reduced, computa- 
tion speed was strongly increased, OpenMP parallelisa- 
tion^ was added where possible, common functions used 
by more than one module were bundled in central libraries, 
and strict adherence to language standards was enforced for 
compiler-independence. 

As a result, there is now a simulation pipeline capable of 
producing realistic data streams for a full year of observation 
with individual Planck detectors within a day of CPU time on 
a computer with sufficiently large main memory («;10GB). For 
near-realistic parameter choices the entire pipeline can success- 
fully be run on standard PCs. Due to the strict adherence to lan- 
guage standards and avoidance of platform-specific libraries, 
the package cleanly compiles on four different types of UNIX 
platforms in both 32- and 64-bit architectures. Furthermore, 
it has pr oven suitable to ru n on a computational grid infra- 
structure haffom et alJ2"on5ll . The pipeline code is freely avail- 
able within the Planck consortium. 

After giving a brief overview of the Planck mission in 
Sect. 12 we describe in Sect. |3 the layout of the pipeline, i.e. 
the decomposition of its simulation tasks into modules, and in 

- ^http://openmp.org^ 
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Sect. I^the data it produces. Sections |5] and |6l list the pipeline 
components common to all modules, and the individual simu- 
lation modules themselves. A summary is presented in Sect.0 

2. Summary of Planck's characteristics 

Here we briefly summarise some characteristics of the Planck 
mission which are of crucial importance for the design and con- 
struction of the simulation pipeline, and for understanding its 
function. 

Planck will observe the microwave sky from a Lissajou or- 
bit centred on the outer (and unstable) Lagrangian point L2 of 
the Sun-Earth/Moon system, 1 .5 x 10^ km from Earth. This en- 
ables the spacecraft to always have the Earth and the Sun at its 
back, and the Moon at a small angle relative to its backward 
direction. Planck will be spin-stabilised, rotating at 1 r.p.m.. 
The optical axis of its Gregorian telescope points at an almost 
right angle with the rotation axis, thus scanning the sky along 
circles. The rotation axis will be kept fixed for 60 minutes and 
then repointed by 2.5 arc minutes such that Planck progresses 
as the Earth moves around the Sun. Planck will thus cover the 
full sky in half a year. The full mission (after reaching L2) is 
expected to last 18 months. 

The rotation is governed by the inertial tensor of the space- 
craft, which is time-dependent due to fuel and coolant con- 
sumption. Repointing the spacecraft every hour causes preces- 
sion which needs to be damped away at the beginning of each 
pointing period. 

The angle between the optical and rotation axes is slighly 
less than 90°. For Planck to cover the full sky, its rotation axis 
will make slow excursions from the Ecliptic plane with an am- 
plitude of several degrees. This scanning pattern implies that 
the sky will be covered inhomogeneously, least densely near 
the Ecliptic and most densely along lozenge-shaped curves 
centred on the Ecliptic poles. Rings on the sky scanned in con- 
secutive pointing periods will intersect, allowing cross calibra- 
tion. The scanning pattern also implies that contact with the 
receiving antenna in Perth (Australia) will be limited to two 
hours per day, during which data down- and uploads will have 
to be completed. 

Two instruments will receive the incoming microwave ra- 
diation. The low-frequency instrument (LFI) has high elec- 
tron mobility transistors (HEMTs) as detectors, allowing mea- 
surement of amplitudes and phases of incoming waves, and 
thus of their intensity and polarisation. It has three channels 
at 30, 44 and 70 GHz. The high-frequency instrument (HFI) 
uses bolometers instead. Four of its channels at 100, 143, 
217 and 353 GHz are polarisation-sensitive, while the highest- 
frequency ones at 545 and 857 GHz are not. Passive cooling 
is insufficient for these instruments. LFI will operate at 20 K, 
while HFI needs to be cooled down to 100 mK. This will be 
achieved by a four-stage cooling system. Residual thermal cou- 
pling between the cooling system, in particular the sorption 
cooler, and the detector electronics, the instruments and the 
telescope gives rise to a noise component which varies with 
the cyclic operation of the cooler. 

The angular resolution of the beams decreases from about 
30 to 5 arc minutes from the 30 to the 857 GHz channels. While 
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Fig. 1. Schematic data flow in a typical Planck simulation 
pipeline. Rectangular components denote parameters or data 
products (see Sect. |4}, whereas elUptic shapes represent mod- 
ules (see Sect.|6}. 



the cores of the beams have approximately Gaussian profiles 
with slightly elliptical cross sections, their wings are extended 
and can receive weak signals from regions of the sky which are 
quite far away from the direction of the beam core. These so- 
called far-side beam lobes thus add a weak noise contribution. 
The feed horns of the individual detectors are slightly inclined 
from their positions in the focal plane towards the optical axis 
of the telescope, which implies that at a given time all detectors 
are observing different points on the sky. 

The channel width is approximately 20% of the central fre- 
quency for LFI, and 30% for HFI. During scans, the signal 
received by the detectors will be integrated over certain time 
intervals and read out at frequencies varying between 30 and 
200 Hz. The sampling frequency needs to be high enough for 
the beam to proceed sufficiently less than its own width along 
the scanning path during one sampling period. 

Slow gain drifts in the detector sensitivity cause Planck's 
sensitivity to vary slowly along the circular scan paths. The fre- 
quency power spectrum of these gain drifts falls approximately 
as a power law in frequency and is thus called 1// noise. This 
noise component gives rise to stripe patterns in the measure- 
ments which follow the scan circles on the sky. 
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3. Pipeline layout 

Fig. [2 illustrates a typical pipeline setup used for simulat- 
ing detector timelines and pointings starting from basic inputs 
like cosmological parameters, microwave emission models for 
galactic and non-galactic components and a scanning strategy. 
A short description of the individual steps is given below: 

- Initially, a CMB power spectrum is generated from a set of 
assumed cosmological parameters, using the CAMB code 
by Lewis and Challinor (see Sect. 16. II Lewis et al. 2000). 

- A particular realisation of this power spectrum is generated 
using the syn_alm_cxx module (Sect. lOt in terms of a set 
of spherical-harmonic coefficients 

- The available set of diffuse foreground component maps is 
transformed into spherical-harmonic coefficients using the 
anafast_cxx module (Sect. l6.2l i. 

- A combined set of «/,„ coefficients is generated using 
the almmixer module, weighting the individual compo- 
nent maps by convolving their electromagnetic spectra 
with the frequency response of the chosen Planck detector 
(Sect.lO. 

- A spherical-harmonic transform of an idealised beam 
template is generated using the gaussbeampol mod- 
ule; alternatively, a realistic, simulated beam pattern can 
be expanded into spherical harmonics using beam2alm 
(Sect.lO. 

- Since a real-space representation of the beam is required 
for efficient convolution with point-sources, the alm2grid 
code (Sect. l6^ is used to generate a beam profile sampled 
on a grid which is equidistant in the polar angles {§, (p). 

- Using the spherical-harmonic coefficients of the combined 
sky and the beam, the totalconvolve_cxx module pro- 
duces a so-called ring set (see Sect. 13, from which the in- 
tensity measured by the detector can be determined for all 
detector pointings and orientations (Sect. (631 . 

- The simmission module produces a time stream of satel- 
lite locations and orientations from input values describ- 
ing the mission orbit parameters, the scanning strategy 
and the dynamical properties of the satellite (Sect. 16. 6> . 
Additionally, simmission can generate data tables con- 
taining the location of the planets as seen from the space- 
craft, which are used by the point-source convolver. 

- The multimod code (Sect. 16. 7> processes the ring sets, 
satellite pointing information and beam profiles, and gen- 
erates time streams containing detector pointings and the 
measured signal. In addition to scanning the convolved dif- 
fuse sky signal (Sect. l6^7?2t . multimod also computes the 
signal of the CMB dipole (Sect. l6.73t . the point-source sig- 
nal (Sect. l6^T4t . the noise generated by the sorption cooler 
on-board the spacecraft (Sect. I6.7.5> . and the instrumen- 
tal 1 //-noise (Sect. I6.7.7> . The sampling and integration 
techniques used for the individual detectors (Sect. I6.7.6> 
are also taken into account. 

Mainly for diagnostic purposes, multimod can also gen- 
erate co-added sky maps (which are produced by sim- 
ply adding every simulated sample to the pixel containing 
the corresponding pointing, and subsequent averaging) and 



coverage maps in HEALPix^^ format, which can in turn be 
converted to Targa^ image files by the map2tga program. 

It is evident that, in this layout, the pipeline becomes in- 
creasingly mission-specific from beginning to end: while the 
calculation of the CMB power spectrum and the creation of a 
set of fl;,„ is completely independent of the mission setup, the 
almmixer module already needs information about the detec- 
tors' frequency response. For the total convolution step, the ge- 
ometrical beam properties must be known as well, and in order 
to produce time-ordered data, a detailed scanning strategy and 
quantities like the satellite's inertial tensor must be specified to 
the simmission module. The multimod program additionally 
requires the detectors' noise characteristics, sampling frequen- 
cies, and integration time constants. 

Conceptually, the simulation process can be subdivided 
into two parts: all modules up to and including the totalcon- 
volver and beam-generation modules produce whole-sky data 
without time dependence; the subsequent modules generate 
time-dependent, localised output. 

It is important to note that the layout described here is by 
no means fixed, but can be adjusted to whatever is needed to 
perform a particular simulation task. This can be achieved by 
using the pipeline editor code, which is part of the Planck pro- 
cess coordinator package. 

The following sections give a more in-depth description of 
the data types and codes constituting the simulation pipeline. 
While this documentation should suffice to decide in which sit- 
uations a given module or data type can be used, it does not 
give detailed usage instructions (such as lists of mandatory pa- 
rameters). This topic is covered by the document "Compilation 
and Usage of the Planck Simulation Modules", available in 
Planck's LiveLink repository or from the authors, which serves 
as a detailed manual for practical use of the simulation pack- 
age. 

4. Data products 

This section documents the various types of data sets produced 
by the Planck simulation package which are intended for fur- 
ther use by the Planck Data Processing Centres. 

Unless otherwise stated, all data products of the Planck 
simulation package containing intensities are stored as antenna 
temperatures measured in Kelvin. Angles are generally stored 
as radians. The simulation pipeline produces maps and pointing 
information in ecliptic coordinates. 

Maps: 

The Planck simulation package uses maps in the HEALPix 
pixelisation scheme to store full-sky information like CMB 
intensity and polarisation, foreground components, sky 
coverage information and co-added maps. In addition to 
single-component maps, three-component maps are used 
to store polarised data in the form of the Stokes I, Q and U 
parameters. 



' http://healpix.jpl.nasa.gov" 
_http://astronomy.swin.edu.au/ pbourke/dataformats/tga/1 
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Spherical-harmonic coefficients (aim): 

A band-Hmited function /(??, tp) with maximum multipole 
order l^.^^ defined on the sphere can be expressed as a set 
of spherical-harmonic coefficients which is implicitly 
defined by 

'max ' 

f{§,<p) = YjYj a,„, f)- (1) 

/=0 m=-l 

These complex coefficients are stored in a three-column 
format, where the first column holds the integer index 
f + I + m + I, which unambiguously encodes both / and m 
because \m\ < I. The second and thkd columns contain the 
real and imaginary parts, respectively, of the coiTespond- 
ing coefficient. Since the Planck simulation package exclu- 
sively works on real-valued maps, coefficients with nega- 
tive m are not stored because of the condition a/ = a^^^ 
for real- valued functions on the sphere. The index column 
is assumed to contain monotonically increasing values, and 
coefficients not listed in the table are assumed to be zero. 
For polarisation information, which can also be stored in 
the form of spherical-harmonic coefficients, the relations 
between real-space and spherical-harmonic space are 
more involved because the Q and U Stokes parameters are 
quantities of spin +2, and must therefore be expanded in 
spin- weigh ted harmonics Their e xact definition can 

be found in lZaldarriaga & SeliakI ( Il997l) . 

Power spectra (C/): 

The power spectra read and generated by the Planck sim- 
ulation modules consist of either one or four components: 
both variants contain the auto-correlation C/j-r of the tem- 
perature fluctuations, the four-component variant addition- 
ally contains the auto-coiTelations Ci^ee and Ci^bb of the E 
and B polarisation modes, as well as the cross-correlation 
CijE between temperature and E-type polarisation. The 
following normalisation conventions are used: 

Ci,aa = ^[^jTJ '"'"''^ (autocorrelation) (2) 

m=-/ 

and 

1 ' 

Ci,ab = (cross correlation) (3) 

m=-/ 

Time-ordered data: 

The time-dependent output signal of the detectors is 
stored as a single column of data. Since data streams 
for the complete mission can become very large, and 
many computer platforms do not support files larger than 
2 GBytes, the Planck simulation modules can optionally 
split time-ordered data into several files. 

Detector pointings: 

This time-ordered quantity is stored in a two- or three- 
column format. The first two columns describe the 



co-latitude and the longitude of the beam centre (in 
radians), while the optional third column contains the 
beam orientation, given as the angle between the tangent 
to the local meridian (pointing northwards) and the x-axis 
of the rotated beam coordinate system. Similarly to the 
time-ordered data, detector pointings can be split into 
several files to avoid the 2-GByte limit of the file size. 

Detector information: 

Several Planck simulation modules need information about 
Planck's detectors such as their position in the focal plane, 
the orientation of their optical axes relative to the tele- 
scope's optical axis, their noise characteristics, spectral re- 
sponses, sampling properties, ellipticities etc. This infor- 
mation is centrally stored and maintained in the so-called 
focal plane database, which was originally a text file, but 
has now been converted to a FITS table for usage within 
the Planck simulation package. Recently, a new tool was 
added which allows extracting the required detector infor- 
mation from HFI's Instrument Model Object (IMO) files 
and converting it to FITS. 

In addition to the detector characteristics, the focal-plane 
database also contains the angle between the nominal spin 
axis of the satellite and the optical axis of the telescope. 

Ring sets: 

This data type is written by the totalconvolver module 
(Sect. 16. St and used by the interpolator module (Sect. 
I6.7.2> . It contains the three-dimensional array 

'„„IX 

r(,^£,0,m")= £ r,„,,„, , (4) 

m,m' 

where 4>k, (t> and Tm„ '.m" are defined in Eq. (8) of 
IWandelt & Gorskil ( 1200 ll) . We do not carry out the Fourier 
transform over the m" direction because the necessary in- 
terpolation can be caiTied out more easily in this represen- 
tation (see Sect. l6.7.2l for details). 

For real-valued skies and beams, the condition 
T{(f)E,<f>,m") = T*{(f)E,<f>,-m") is fulfilled, thus the 
coefficients with negative m" and the imaginary part of 
T{(pE, 0, 0) need not be stored. 

5. Common components 

5.1. External libraries 

5.1.1. FITS I/O 

Since all mod ules are presently cor nmunicating via FITS files 
(lNOSlll999l) . the cf itsio Ubrary jPencel 19991) is required for 
the pipeline to operate. Currently, version 2.5 10 of the library is 
used. In the future, it is planned to store the data in a (possibly 
object-oriented) data-base system. 

5.1.2. FFT 

Several modules need to perform fast Fourier transformations, 
which can be divided into two diff"erent usage scenarios: 
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- Most codes (like the transformations of maps between 
real and spherical-harmonic space carried out within 
the HEALPix package) perform FFTs of many different 
lengths, but each FFT of a given length is only performed 
a few times. Here the Planck simulation package contains 
a port to the C language of the FF TPACK code originally 
developed bv IS warztraubeJ ( Il982t) . If the length of the re- 
quired F FT has very l arge prime factors, Bluestein's algo- 
rithm teluestein 1968) is employed; the choice of the algo- 
rithm occurs automatically. 

Both complex- and real-valued transforms are supported. 

- Several other, still experimental, modules require a highly 
optimised FFT transform of a given length, which is exe- 
cuted many times. This task is handled by the FFTW library 
jprigo & .Tohns on 1998), of which the Planck simulation 
package currently uses version 3.0.1. 

5.2. Other common functions 

Apart from using externally developed general-purpose soft- 
ware such as libcfitsio, the codes have many additional 
commonalities: they perform file input and output in a very 
similar way (using only a small subset of the cfitsio func- 
tion), they deal with a common set of data types (see Sect. |^, 
and they need to operate on them in similar ways. Keeping 
a separate implementation of these facilities for each module 
would cause unnecessary code duplication, leading to software 
that is both eiTor-prone and hard to maintain. For this reason, 
common functions are implemented exactly once and linked to 
all programs that need it. However, since some modules in the 
Planck simulation package are written in Fortran 90 and others 
in C++, a few common components (most notably the support 
functions concerned with data I/O) had to be implemented in 
both languages to allow convenient usage. 

5.2.1. Random number generator 

Several modules, most notably the 1 //-noise generator and 
the syn_alin_cxx module, which creates sets of spherical- 
harmonic coefficients from a power spectrum, require streams 
of random numbers. Using the default random-number genera- 
tors supplied with the programming language is not appropriate 
in our situation, since the period length and overall quality of 
these random numbers depend on the operating system and the 
compiler, whereas the data products created with the Planck 
simulation package should be identical on all supported plat- 
forms (for identical random seeds). 

Currently, the P lanck simulation modules use a generator 
of the xorshift type (lMarsaglial2003l) with a period of 2'^^ - 1, 
which is more than sufficient for the necessary calculations. In 
addition to its excellent randomness properties, the algorithm 
is also very efficient. 

Since most codes require random numbers with a Gaussian 
distribution instead of uniform ones, the simulation pack- 
age contains an imp lementation of the Box-Muller method 
ilBox&MulleHll958h for generating the former from the lat- 
ter. 



5.2.2. Input of simulation parameters 

All Planck simulation modules follow a single convention for 
how simulation parameters (like cosmological parameters, map 
resolution, detector names, names of input and output files 
etc.) are obtained from the environment. This is done via plain 
ASCII files obeying a simple format. Each parameter is defined 
in a separate line using the form: 

[name] = [value] 

Empty lines are skipped. Lines starting with a hash (#) sign are 
ignored and can thus be used as comment lines. 

5.2.3. Handling of data I/O 

Data exchange between the various pipeline modules is per- 
formed exclusively via FITS files. Although the cfitsio li- 
brary provides a comprehensive interface for creating and read- 
ing these files, calling it directly from the modules is in most 
cases more complicated and error-prone than necessary for the 
limited set of data types used within the Planck simulation 
package. To simplify data input and output for the module au- 
thor, a powerful but rather easy-to-use set of "wrapper" func- 
tions around cfitsio was written in both Fortran 90 and C++, 
which supports all input and output operations performed by 
the Planck simulation modules. 

This abstraction from the interface of the underlying library 
has the additional advantage that the external data format can 
be changed with relatively little impact on module code, simply 
by adjusting the implementation (without modifying the user- 
visible interface) of the common input-output functions men- 
tioned above to the new format. The anticipated switch from 
FITS files to the Planck Data Management Component as a 
means for transporting data between modules will greatly ben- 
efit from this design. 

On an even higher level, functions exist which read and 
write entire maps, spherical-harmonic coefficient sets and 
power spectrum objects; this code is also used by many mod- 
ules. 



6. Module descriptions 

In the following subsections, we describe the function of all 
modules in the Planck simulation package. Where possible, we 
only give a short description of the underlying algorithms and 
cite the original works introducing the numerical methods; if no 
such publication exists, we provide a more detailed description. 
For the most resource-consuming applications we also spec- 
ify their CPU, memory and disk-space requirements as well 
as their scaling behaviour with regard to relevant input param- 
eters (such as the map resolution or the maximum multipole 
order Z,nax)- Parallelisation of a code is also mentioned where 
apphcable. 
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6.1. CAMB 

For generating CMB power spectra, the November 2004 ver- 
sion of the CAMB^ code (Lewis et al. 2000) is used. This code 
is based on CMBfast ( Sehak & Zaldarriag a 1996), but is more 
portable and can be more easily integrated into the Planck sim- 
ulation package, mainly because it does not rely on unformat- 
ted files for storage of intermediate data products. 

Memory and CPU requirements of the CAMB code 
are rather modest for the parameter sets typically used 
for Planck simulations (memory consumption < 100MB, run- 
time < 1 minute), hence it can be run without any problems on 
desktop machines. 

6.2. HEALPix 

The Planck simulation pipeline m akes extensive use of the 
HEALPix pixelisation of the sphere dOorski et alJ2002lE005i) . 
It is implemented in a software package which contains many 
associated algorithms, among them some exploiting the effi- 
ciency of the HEALPix scheme for spherical-harmonic trans- 
forms. The main usage areas of the package in the simulation 
pipeline are: 

- generating spherical-harmonic coefficients from power 
spectra Q; 

- conversion between aim and full-sky maps; 

- generation of co-added maps from time-ordered data; 

- generation of efficient look-up tables for the point-source 
convolver 

A large part of the HEALPix software package, which is 
available in Fortran 90, was re-implemented in C-I--1- in order 
to be conveniently accessible from within the C++ simulation 
modules. During this migration, the function of the synfast 
tool, which essentially generates sky maps from power spectra, 
was split into two separate facilities, syn_alin_cxx (which gen- 
erates aim from a power spectrum) and alm2map_cxx (which 
transforms o/,,, coefficients to a HEALPix map). As can be seen 
in Fig. n this allows feeding the results of syn_alin_cxx di- 
rectly into the alrmnixer module (Sect. 16. 3> . and avoids the 
conversion to real space and back to «/„,, which was neces- 
sary in earlier versions of the Planck simulation package and 
caused some avoidable loss of accuracy and waste of computer 
resources. 

The ported code (especially the conversions between maps 
and fl/,„ coefficients) was optimised with regard to memory 
and CPU time consumption. On shared-memory computers, 
OpenMP parallelisation is used. 

An addition to the HEALPix facilities is the alm2grid 
program. This code converts «/,„ coefficients to values on an 
equidistant grid in the polar angles {&,(p), whose extent in & 
can be chosen freely. Using this format, it is easy to calculate 
a high-resolution real-space representation of a beam pattern, 
which in turn is a necessary input for the point-source con- 
volver 

The C++ package also contains a testsuite to validate the 
correctness and accuracy of the ported code. 

' |http://camb.info. 



Starting with version 2.0, the official HEALPix package 
contains the described C++ modules, and is distributed under 
the GNU General Public License^. 

6.3. skymixer/almmixer 

The skymixer component serves two purposes: 

- It converts maps of the CMB and diffuse foreground 
sources (such as the Galactic synchrotron and dust emission 
and the thermal and kinetic Sunyaev-Zel'dovich effects) 
from the units in which they are typically provided (e.g. 
MJy/sr at a certain frequency for the Galactic maps, and the 
Compton y-parameter for the thermal Sunyaev-Zel'dovich 
maps) to antenna temperature observed by a given Planck 
detector During this conversion, the frequency response of 
Planck's detectors is taken into account, i.e. the frequency- 
dependent sky signal is convolved with the sensitivity of the 
detectors not only at the nominal detector frequency, but in 
a frequency band around it. Currently, the spectral detector 
response can have Gaussian, top hat or 6 form. 

- It adds the resulting antenna-temperature maps to a single 
map with all selected components. 

The currently employed code is a re-implementation of a 
version by R. Ansari, which was based on the SOPHYA^ li- 
brary. 

The microwave sources currently supported by skymixer 

are: 

- CMB given as temperature fluctuations in Kelvin: 

This kind of map is typically generated from a power spec- 
trum Ci produced by the CAMB code. 

- Galactic synchrotron emission, given in MJy/sr at 
408 MHz: 

The template was provided by G. Giardino; de t ails co ncern- 
ing its creation are given by 'Giardino et alJ ( l2002h . This 
component is extrapolated to the Planck frequency bands 
by using a power-law spectrum with an exponent of -0.75 
between 408 MHz and 22 GHz, and with an exponent of 
-1.25 above 22 GHz. 

A spectral-index map derived from Wilkinson-MAP data is 
also available, but has not yet been included into the extrap- 
olation algorithm. 

- Galactic dust emission, given in MJy/sr at 100/im: 
Currently this source is handled by a one-component ap- 
proximation provided by C. B accigalupi. The full tw o- 
component model presented by iFinkbeiner et alJ (Il999l) is 
currently being implemented. 

The employed te mplate map was published by 



ISchlegel etalJ(IT99^ 

- Galactic free-free emission: 

This map was derived from a template rn ap of galactic Hq.- 
emission published bvlFinkbeineJ ( |2003l , according to the 
model of lVaUs-Gabaudl(ll998l) . 
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- thermal and kinetic SZ-efFect, given as the Compton-y and 
Compton-w parameter, respectively: 

Briefly, these maps were produced using the galaxy-cluster 
distribution in space, peculiar velocity and redshift found 
in the Hubble Vol ume simulation {Jenkins et al. 2001; 
IColberg etal] l2000l) . attaching appropriately scaled gas- 
dynamical cluster simulations to the cluster positions and 
projecting the Compton-y and w parameters on the sphere. 
Details about the template maps are given bv lSchafer et alJ 
l2004a). 

Additional maps for the thermal and kinetic SZ-efFect of 
the local super cluster structure have become available re- 
cently. In this case, the full sky SZ maps are directly cal- 
culated from a hyd odynamical, con strained simulation of 
the local universe (Dolag et al. 2005), which spans a sphere 
of 1 lOMpc radius. These maps are complemented by con- 
sidering the galaxy-cluster distribution outside this sphere 
using the Hubble Volume simulation as described before. 

- Emission by rotating CO molecules: 

This foreground component only contributes to the HFI 
channels. The template map was published by Dame et al. 
1 1996, 200 JJ. We use the method given bv .Schafer et al.. 
1 2004b) to determine the intensities in Planck's frequency 
bands. 

Since the output of the skymixer is most frequently used 
in a subsequent run of the totalconvolve_cxx module (see 
Sect, \6.5i which requires a/,,, coefficients as input, and since 
the synthesised CMB signal can be very efliciently created 
in the form of spherical-harmonic coeflicients o/,,,, an addi- 
tional variant of the skymixer code was implemented. This 
code (named almmixer) performs exactly the same task as 
skymixer, but operates entirely in spherical-harmonic space. 
Using this code eliminates the need for converting between an- 
gular and spherical-harmonic space twice during a run of a typ- 
ical simulation pipeline and thus saves a considerable amount 
of CPU time. 

Both codes allow the suppression of the CMB monopole 
in the output, which is usually not needed by any subsequent 
analysis, and which would decrease the numerical accuracy of 
the output. 

The runtime of both programs is rather short and dominated 
by the input and output operations. The memory consumption 
of the skymixer is low (below 10 MBytes) and independent 
of the map size, since it processes them chunk by chunk. In 
contrast, the almmixer module holds two o/,,, sets in memory, 
so that it requires approximately 16/^^^ bytes of storage. 

6.4. Beam generation 

Several Planck simulation modules require information about 
the beam shapes of Planck's detectors, typically in the form of 
a set of spherical-harmonic coefficients. 

One way to obtain these is the module gaussbeampol, 
which takes a detector name as input, looks up the detector's 
beam characteristics in the focal-plane database, and produces 
the corresponding coefficients. The gaussbeampol module 
can generate axisymmetric or elliptical Gaussian beams with 



and without polarisation, and allows arbitrary orientations for 
polarisation and ellipticity. 

Realistic beam patterns are often delivered in the GRASPS 
format**, which in turn consists of the GRD and CUT sub- 
formats. The beam2alm code can convert beam files of both 
sub-formats to ai,,, coefficients, and also allows combining two 
GRASPS files (containing the full and main beam, respec- 
tively) to a single set. 

The point-source convolver requires a real-space represen- 
tation of the beam, which can be obtained from the ai„ coeffi- 
cients using the alm2gridor alm2map_cxx modules. 

6.5. Total convolution 

This module takes as input the spherical-harmonic coefficients 
of a sky map and a beam, and computes the convolution of 
sky and beam for all possible directions (i?, tp) and orientations 
(tfr) of the beam relative to the sky. Both unpolarised and po- 
larised convolutions are supported. In the polarised case, three 
sets of spherical harmonics are required for sky and beam, re- 
spectively, for their Stokes-I, Q, and U parameters. The output 
then consists of the sum of the three separate convolutions. 

Our i mpleme ntation is analogous to the algorithm pre- 
sented bv IWandelt & Go rski (200 it, inc l udin g the extensions 
for polarisation given by Challino r et alJ (EqOO). However, the 
code was modified to allow complete shared-memory paralleli- 
sation of all CPU-intensive tasks (starting from the partial par- 
allelisation contributed by Stephane Colombi), and the sym- 
metri es inherent in the matrix d'^^^, (introduced in Eq. (6) of 
(Wandelt & G6rski200ll) were used for further CPU and mem- 
ory savings. With the new code, convolutions up to /^ax = 1500 
and m,nax = 2 can be done in less than five minutes on a modern 
desktop computer 

The resource consumption of the totalconvolve_cxx 
module scales roughly with /max'"max for memory and 
^max '"max for CPU time. Although the code is fully parallelised, 
it will most likely not scale very well, since the cache re-use of 
the algorithm is poor and a saturation of the memory bandwidth 
will limit overall performance. 

The output consists of the three-dimensional complex- 
valued array T{i9, ip, m), which contains Zmax + 1 + x points in 
)?-direction (x is a small number of additional points needed 
for interpolation purposes), 2Zn,ax + 1 points in (^-direction, and 
'«max + 1 values in m-direction. In contrast to the original im- 
plementation by B. Wandelt, the final EFT is only carried out 
over &- and (^-directions, because interpolation is easier for this 
case; see also Sect. 16.7.21 

To increase the spatial resolution of the output ring set, the 
working array can be zero-padded to any desired /max.out ^ ^max 
before the EET is performed. 

6.6. Satellite dynamics 

The Planck simulation package makes use of the simmission 
code developed by Eloor van Leeuwen and Daniel Mortlock 
lyan Leeuwen et al. 2002: Challinor et aL..2002) . which takes 

^ ihttp://www.ticra.com| 
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as input a wide variety of mission-related parameters like the 
exact orbit around the outer Lagrangian point L2 of the Sun- 
Earth/Moon system, the scanning strategy, the satellite's iner- 
tial tensor etc. 

The most important output of the siitimission module is 
a table containing the position and orientation of the satellite 
in fixed time intervals during the entire mission. In addition, it 
optionally calculates the positions of the Sun and planets rela- 
tive to the Planck spacecraft, which is important for the point- 
source convolver. 

6.7. multimod 

The multimod code comprises most of the functions offered 
by the Planck simulation package for creating and manipulat- 
ing time-ordered information. This includes creation of detec- 
tor pointings and signal time streams from various sources like 
the CMB dipole, point sources and detector noise. 

At first glance, this approach appears to contradict the goal 
of modularity stated at the beginning of this paper, which would 
suggest writing many small pieces of code implementing a 
single task each. This strategy, however, would have a nega- 
tive impact on performance, since in most of the simulation 
codes working on time-ordered data, the disk input or output 
times are at least comparable to the computations themselves. 
Combining all the effects in one program allows to read and 
write the data only once and perform all manipulations without 
storing intermediate data products on disk. 

Even though all the functions are linked together into one 
executable, the various components are still implemented in 
a very modular fashion, so that they can be easily used to 
build smaller modules with reduced generality. For example, 
the Planck simulation package contains a stand-alone point- 
source convolution code, which includes exactly the same C++ 
classes also used by the multimod code. 

The multimod program processes data in chunks of one 
pointing period at a time, which corresponds to an hour of mea- 
surement. This keeps the amount of main memory occupied by 
time-ordered data rather low. On the other hand, a pointing pe- 
riod still contains enough individual samples to allow efficient 
OpenMP-parallelisation of the code. 

Schematically, multimod performs the following steps for 
each pointing period: 

1. Calculate the detector pointings at a sampling rate of 
ringres/60s from the satellite pointings (see Sect. 16. 7.11 . 

2. Calculate and add all of the requested signal contributions 
(except for l//-noise) at a sampling rate of ringres/60s 
(see Sections l6T2l - l6.7.5> . 

3. Process the calculated signal in the sampler module (see 
Sect. l6T6l l. 

4. Add the 1 //-noise if requested (see Sect. l6.7.7l i. 

5. Produce the detector pointings at the detector sampling 
rate. 

6. Write the requested signal and pointing streams to disk. 

7. Update the coverage and coadded sky maps. 

Each of these steps is only executed if its output was requested 
by the user or if it is required by subsequent steps; i.e. if the 



user only asks for the generation of detector pointings, only 
steps|5land (partially)|6lare executed. 

The parameter ringres can be chosen freely by the user; 
it determines the sampling rate used for producing the "ideal" 
detector signal (i.e. the signal before the modification by the 
detector electronics and the sampling process). In theory this 
sampling rate should be infinitely high; for practical purposes 
it should be sufficient to choose ringres/60s x /s^mp for 
smooth signals (like the CMB and galactic foregrounds) and 
ringres/60s ~3/samp for the SZ and point-source signals. 

6.7.1 . Detector pointings 

Using the time stream M„ of satellite orientations generated 
by the simmission module, the detector-pointing component 
produces time-ordered detector pointings at arbitrary sampling 
rates. For this purpose, it is necessary to determine the satel- 
lite orientation at any given point in time. This is achieved by 
extracting the rotation matrix R„ between every M„ and M„+i. 
From this matrix, an axis r„ and angle a„ of rotation can be 
computed. The satellite orientation at a time t with t„ < t < t„+i 
can then be approximated by 

M, = RM (r„, a„] M„ , (5) 

\ t„+i - 1„ j 

where RM(r, a) is a rotation matrix performing a rotation of a 
around the axis r. 

After this rotation interpolation another rotation matrix 
must be applied, which describes the detector position and ori- 
entation relative to the satellite coordinate frame. This matrix 
depends on the angle between the optical axis and the nominal 
spin axis, as well as the detector-specific angles (^,„,, j?,,,. and t/r,„, 
given in the focal-plane database, which specify the position 
and orientation of all detector horns relative to the coordinate 
frame of the focal plane. 

6.7.2. Interpolator 

The interpolator module's task is to determine the radiation 
intensity falling onto a detector for a given co-latitude d-, longi- 
tude (f and beam orientation t/^. To accomplish this, the output 
of the totalconvolve_cxx modul e for the desired sky signal 
and the detector's beam is used. As IWandelt & Gorskj ( l200ll) 
point out, this dataset is sufficient (for a given /max and mmax 
of the beam) to reconstruct the exact signal for any (??, ip, xf/)- 
triple. However, exact evaluation is prohibitively expensive in 
practice, since it would involve a sum over Z^ax '"max terms for 
a single data point. 

The approach used in the interpolator is based on the 
assumption that /^ax will be rather large 4000) for realistic 
simulations, while the beam's ;«max is usually not larger than 20 
(i.e. the beam pattern has no high-frequency azimuthal asym- 
metries). The high /,nax implies a fairly high resolution in and 
for the totalconvolver's output, so that polynomial interpola- 
tion (whose order can be chosen via a parameter) can be used 
for those two directions. In the (/^-direction, this approach is 
not accurate enough because of the sparse sampling, so that the 
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Fourier sum of the i/r-components must be calculated: 

'"max 

2 r,„terpol(>?,^,m)e™^ (6) 

A straightforward calculation of e""'^ would be quite CPU- 
intensive because of the required trigonometric function calls, 
but making use of the trivial recursion relation e'('"+i)'/' = 
gimiA giiA reduces calculation time significantly. This optimisa- 
tion is only applicable here because m,nax cannot become very 
large; if it did, the roundoff errors accumulated in the recursion 
could destroy the accuracy of the result. 

The memory consumption of the module is proportional to 
'max '"max, its runtime scales with m^ax "samples- Since interpo- 
lations for different pointings are independent, parallelisation 
was trivial. 

By default, the interpolator module calculates the total 
intensity observed by the detector Optionally, any of the indi- 
vidual Stokes parameters 1, Q and U can also be extracted, but 
this calculation is only exact for axisymmetric polarised beams. 

When using linear interpolation, it has been observed that 
the power spectrum of the resulting co-added maps was sup- 
pressed in comparison to the input maps at high I. This is a 
consequence of the finite resolution of the ring set, combined 
with the non-perfect interpolation in i? and (p. Although this 
effect cannot be eliminated entirely, increasing the interpola- 
tion order and/or increasing Zmax.out in the totalconcolver 
reduces the problem dramatically. 

6.7.3. Dipole 

This component calculates the time-dependent dipole signal 
observed by a given detector The user can choose whether 
the kinematic dipole component (caused by the motion of the 
satellite in the Solar System, with the Solar System within the 
Galaxy, with the Galaxy within the Local Group, and so forth) 
should be included, whether relativistic effects should be calcu- 
lated, whether only higher order corrections should be returned 
etc. The default output quantity is the dipole signal in antenna 
temperature, but thermodynamic temperature can also be cal- 
culated. 

6.7.4. Point-source convolver 

The task of the point source convolver is to simulate time 
streams of the signal created by fixed and moving point 
sources. For fixed sources, this could in principle be achieved 
also by adding another foreground to skymixer or almmixer, 
but this is not very practical because the maximum multipole 
moment required in the totalconvoler to accurately treat 
this kind of strongly localised sources would be very high. For 
moving sources, this solution is not viable at all, since all sig- 
nal contributions entering the mixing modules must be constant 
over the mission time scale. 

For these reasons, the point source convolver operates in 
real space, taking as input the detector pointings and convolv- 
ing the sources close to a given pointing with the properly ori- 
ented beam pattern. Obviously, an efficient retrieval of sources 



in the vicintity of the beam is crucial for the performance of 
the module. This was achieved by presorting the point source 
catalog into a HEALPix map with low resulution (A^si(je=16), 
and by identifying, for all detector pointings, the nearby point 
sources using the query_disc() function of HEALPix with a 
radius of several times the beam FWHM. 

While fixed point sources stay in this HEALPix lookup ta- 
ble during a whole simulation run, the positions of the moving 
point sources are calculated once per pointing period. These 
positions are then sorted into the lookup table accordingly and 
taken out again after the signal for the pointing period has been 
calculated. 

As mentioned above, it is recommended to choose a high 
value for the ringres parameter when calculating the point 
source signal, because of their very localised nature. 

For the fixed point sources, the radiation flux in each 
Planck frequency band is taken from the point source cat- 
alog. For moving point sources (i.e. planets and asteroids), 
the corresponding positions in time are calculated using the 
JPL HORIZONS system^, whereas the radiation intensity of 
the moving point sources is calculated usin g an extended 
Wrig ht & Odenwald th ermal emission model ( Wright '197^; 
N eu^ebauer et alJl97ll) which has been adapted separately for 
rocky planets, gaseous planets and asteroids, and which ac- 
counts for effects like distance to the satellite, illumination ge- 
ometry, beaming, surface roughness, heating and cooling of the 
surface, etc. This will be described in more det ail in a sepa- 
rate publication (R. Hell, in preparation; see also lSchafer et alJ 
l2n04hi. 

6.7.5. Sorption-cooler noise 

Temperature fluctuations in the detectors caused by the oper- 
ation of Planck's sorption cooler constitute an important sys- 
tematic effect and have to be included in realistic simulations 
of time-ordered data. To this end, the Planck simulation pack- 
age contains the glissando code contributed by LFl, which 
models this effect and will be fully integrated with multimod 
in the future. Since the LFI and HFl detectors may react quite 
differently to temperature fluctuations, it is likely that a com- 
pletely different code will ultimately be required for the HFI 
channels. 

6.7.6. Sampler 

The sampler component takes as input the idealised signal im- 
pinging on the detector, which is produced by the components 
described above, and applies to it a model of the detector's mea- 
surement process. Two effects play an important role: 

- Every sample recorded by any detector is in reality the av- 
erage of several /flif samples taken in the time since the last 
slow sample. 

- All HFl detectors have an exponential impulse response 
with a time constant Tboi, so that each fast sample must in 

' Jittp://ssd.jpl.nasa.gov/horizons.htnif| 
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turn be calculated via the integral 



«bol(0 = I t/f'isky(f') 



exp[ -(f - f')/Tboi] 



(7) 



together with a process with tq 
white noise no(0- 

The resulting SDE spectra 



for the high-frequency 



The algorithm approximates this integral by taking un- 
equally spaced samples: 



Pi(f)^{Tf+{2 7Tffy 



(11) 



Shol(t) ~ ^ «sky(4) ' 



(8) 



can be combined to an approximation of the target spectrum 
P*{f) (Eq.|5} by combining the individual processes with prop- 
erly chosen weights c,-: 



where t', - t + Tboi ln(^/A^). 



n(t) = ^ C; Xi(t) , 



A more detailed descrip tion of the implementation can be P„(f) - 'S^ Piif) ~ P*{f)- 
found in lGrivell & Manrj [l999). " ; ' ' 



(12) 



(13) 



It should be mentioned that both effects introduce a time 
delay into the sampled time-ordered data, i.e. that the signal 
written with the time stamp f„ was produced entirely from mea- 
surements at f < t„, which in the case of HFl reflect a signal 
that was seen even earlier, because of the bolometer's time re- 
sponse. This must be taken into account when time-ordered sig- 
nal data and detector pointings are used to produce sky maps. 

The sampler code can be trivially parallelised and has per- 
fect cache behaviour, so that it should scale very well on mul- 
tiprocessor machines. 

As expected, the sampler module produces data at a rate 
equal to the detector's sampling frequency. 

6.7.7. l//-noise 

The electronics of both LFI and HFI detectors produces noise 
which has a power-law spectrum with spectral index < y < 2 
between a lower frequency and a higher (or knee) fre- 
quency/knee- Below and above those frequencies, the spectrum 
is white. The targeted spectral power is thus 



p:(/) = p„(i + (///™„)V + p. 



white 



(9) 



where f white is the high-frequency, white-noise level, and Pq = 
-Pwhite (/kneeZ/min)'' - 1)- The valuc of y is » 1.7 for LFI and 2 
for HFl detectors. 

A noise time-stream with this spectrum could be obtained 
by modelling it in Fourier space with random coefficients and 
then performing an FFT. However, considering the fact that the 
coherence length for an HFl detector would be approximately 
10** samples, this method is rather expensive in terms of main 
memory consumption. 

The noise is more efficiently generated by sampling and 
adding numerical solutions of stochastic differential equations 
(SDEs) of die form 



Xi(t) = -Xi(t)/Ti + ^i(t) 



(10) 



where r,- is the autocorrelation time of process x,(f). The ^,(f) 
are individual white noise processes which are implemented in 
the discretised form of Eq. (llOt as normalised series of inde- 
pendent Gaussian random numbers ^,(f) with autocorrelation 
function <^,(f)^/(i)) = 6(t - s). 

In order to have a good logarithmic frequency coverage, 
a logarithmic grid of r between 2;r//knee and 2n/ f„^i„ is used 



More technical detail on choosing the weights, initialising the 
SDE processes, and optimising the computational performance 
will be covered in a separate publication (EnBlin et al., in prepa- 
ration), but see alsolKeihanen et al. (2004). 

Since the output of the noise generator is added to the out- 
put of the sampler module, both are naturally produced at the 
same sampling frequency. 

Parallelisation of the algorithm is not straightforward, be- 
cause every noise sample can only be computed after all of 
its predecessors have been calculated. This prohibits the most 
straightforward strategy of dividing the requested samples be- 
tween multiple processors. In our approach, each SDE process 
calculates its associated time stream separately (which can be 
done in parallel processes), and all of these time streams are 
added at the end. While this solution can be significantly faster 
than a scalar implementation, its scalability is limited by the 
number of SDE processes («; 10). Increasing the number of pro- 
cessors beyond this value achieves no further speedup. 

7. Summary and outlook 

We have described in this paper the current state of the Planck 
simulation pipeline. This package allows a complete simulation 
of the time-ordered data streams produced by the satellite's de- 
tectors, taking as input a set of cosmological parameters, var- 
ious foreground emission models, a satellite scanning strategy 
and parameters characterising the detectors and on-board elec- 
tronics. The availability of such data sets allows extensive test- 
ing of the data analysis software for the mission already be- 
fore launch, which greatly decreases the time interval between 
the availability of raw data and the publication of first results. 
Furthermore the simulation modules can be conveniently used 
to study the influence of assumed systematic efifects or changes 
in the mission parameters on the quality of the final data. 

Obviously, many of the package's modules are stiU under 
development, and additions of new modules are expected. The 
integration with the Planck Data Management Component and 
the mission's Integrated Data and Information System (IDIS) 
will be another important milestone. Any significant changes 
and enhancements will be presented in a follow-up paper. 

The non-f ZancA;-specific modules will hopefully be of use 
outside the Planck collaboration as well. This covers almost 
all of the presented codes, with the exception of the modules 
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for the satellite dynamics and the sorption-cooler noise. A fu- 
ture public release of the non-P/ancfc-specific parts of the pack- 
age is envisaged by the authors. Currently, a partial simulation 
pipeline, which is hosted by the German Astrophysical Virtual 
Observatory, can be tested interactively over the Internet'". 
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